AD-A212  839 


Nava!  Research  Laboratory 

Washington,  DC  20375-5000 


NRL  Memorandum  Report  6478 


Numerical  Evaluation  of  the  Far  Field  Wave 
Pattern  of  the  Radiation  Green’s  Function 

Henry  T.  Wang*  and  E.  Wade  Miner** 

* Center  for  Fluid/Structure  Interactions 
** Center  for  Fluid/Dynamic  Developments 
Laboratory  for  Computational  Physics  and  Fluid  Dynamics 


September  5,  1989 


Approved  for  public  release;  distribution  unlimited 

89  9  26  162 


SECURITY  ClASS'F'CATlQN  QC  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


form  Approved 
OMB  No  0704  0188 


la  REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


2a  SECURITY  CLASSIFICATION  AUTHORITY 


2b  DECLASSIFICATION /DOWNGRADING  SCHEDULE 


4  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

NRL  Memorandum  Report  6478 


lb  RESTRICTIVE  MARKINGS 


3  DISTRIBUTION /AVAILABILITY  OF  REPORT 

Approved  for  public  release;  distribution 
unlimited . 


5  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


6a  NAME  OF  PERFORMING  ORGANIZATION 

Naval  Research  Laboratory 


6c  ADDRESS  (City,  State,  and  ZIP  Code) 

Washington,  DC  20375-5000 


8a  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION 

Office  of  Naval  Research 


8c  ADDRESS  (City,  State,  and  ZIP  Code) 

Arlington,  VA  22217 


6b  OFFICE  SYMBOL  7a  NAME  OF  MONITORING  ORGANIZATION 
(if  applicable) 


Code  4420 


7b  ADDRESS  (City.  State,  and  ZIP  Code) 


8b  OFFICE  SYMBOL  9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
(If  applicable) 


10  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 
ELEMENT  NO 

61153N 


PROJECT 

NO 


WORK  UNIT 
ACCESSION  NO 

DN158-17 


1 1  TITLE  (Include  Security  Classification) 

Numerical  Evaluation  of  the  Far  Field  Wave  Pattern  of  the  Radiation  Green's  Function 


12  PERSONAL  AUTHOR(S) 

Wang,  H.T.  and  Miner,  E.W. 


13a  TYPE  OF  REPORT  13b  TIME  COVERED  114  DATE  OF  REPORT  {Year.  Month,  Day) 

Interim  from  1/87  to  10/88  I  1989  September  5 


16  SUPPLEMENTARY  NOTATION 


COSATI  CODES 


SUB-GROUP 


18  SUBJECT  TERMS  l Continue  on  reverse  if  necessary  and  identify  by  block  number) 

Far  field  wave  pattern  Translating  oscillating  source 
Radiation  Green's  function 
Computer  program  RADGREEN 


19  ABSTRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

Computer  program  RADGREEN  has  been  developed  to  calculate  the  far  field  wave  pattern  of  the  radia¬ 
tion  Green’s  function  corresponding  to  an  oscillating  source  translating  below  the  free  surface.  The  program 
has  the  capability  of  calculating  the  Green’s  function  and  its  derivatives,  or  the  associated  wave  elevation  and 
slopes.  The  far  field  analysis  is  based  on  an  asymptotic  formulation  which  neglects  terms  which  are  of  order 
1/R  and  higher.  Techniques  for  removing  and  attenuating  singularities  in  the  formulation  and  eliminating 
short  noise-like  waves  are  described  in  some  detail.  A  number  of  three-dimensional  and  contour  plots  are 
given  to  show  sample  radiation  wave  patterns  for  a  single  source  and  for  a  line  of  sources  whose  strengths  are 
adjusted  to  approximate  those  for  a  ship  hull  undergoing  heave  or  pitch  motions. 


20  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT 
C3 UNCLASSIFIED/UNLIMITED  □  SAME  AS  RPT 


22a  NAME  OF  RESPONSIBLE  INDIVIDUAL 

Henrv  T.  Wan 


OD  Form  1473,  JUN  86 


21  ABSTRACT  SECURITY  CLASSIFICATION 

■clSers  UNCLASSIFIED 


22b  TELEPHONE  (Include  Area  Code)  22 c  OFFICE  SyMBO. 

202)767-2516  Code  4420 


Previous  editions  are  obsolete 

S/N  0 102-LF-0 14-6603 


SECURitv  QASS'fiC a;  O tv  Ol  this 


CONTENTS 


1.  INTRODUCTION  .  1 

2.  DESCRIPTION  OF  FAR  FIELD  FORMULATION  .  2 

3.  COMPUTER  IMPLEMENTATION  OF  THEORETICAL  FORMULATION  .  6 

4.  NUMERICAL  RESULTS  .  10 

5.  SUMMARY  .  11 

6.  ACKNOWLEDGMENT .  12 

7.  REFERENCES  . .  12 


/w 


NUMERICAL  EVALUATION  OF  THE  FAR  FIELD  WAVE  PATTERN 
OF  THE  RADIATION  GREEN’S  FUNCTION 


1.  INTRODUCTION 

The  method  of  Green’s  functions  is  widely  used  to  analyze  the  flow  around  surface  ship  hulls. 
If  one  neglects  very  short  waves,  such  that  surface  tension  and  viscosity  effects  may  be  neglected, 
these  functions  represent  potentials  of  point  sources  which  are  chosen  to  satisfy  all  the  required  condi¬ 
tions  of  the  problem  with  the  exception  of  the  kinematic  condition  of  no  flow  through  the  ship  hull. 
This  latter  condition,  in  turn  is  used  to  determine  the  strengths  of  the  sources  which  are  placed  on  the 
hull  surface. 

A  critical  step  in  this  procedure  is  the  efficient  evaluation  of  the  near  field  behavior  of  the 
Green’s  function  to  determine  the  mutual  influence  of  the  sources.  Once  these  influences  have  been 
determined,  the  problem  reduces  to  a  system  of  simultaneous  linear  equations  for  the  source  strengths. 
The  near  field  behavior  is  also  of  interest  in  obtaining  the  pressure  distribution  on  the  hull  surface,  as 
well  as  the  bow  wave  in  such  applications  as  slamming  and  deck  wetness. 

The  far  field  analysis  is  often  used  to  obtain  the  damping  and  resistance  forces  due  to  the  waves 
generated  by  the  motions  of  the  ship.  This  usually  leads  to  a  more  accurate  determination  than  the 
direct  integration  of  the  hull  pressure  distribution  due  to  the  greater  simplicity  and  higher  resolution  of 
the  far  field  analysis.  More  recently,  with  the  advent  of  more  accurate  remote  sensing  techniques,  the 
far  field  wave  pattern  is  also  of  interest  from  the  detection  point  of  view. 

The  three  most  widely  used  Green’s  functions,  listed  in  increasing  order  of  complexity,  are  the 
ring  wave  (or  zero  speed)  function  for  an  oscillating  stationary  source,  the  Kelvin  (or  wave  resistance) 
function  for  a  nonoscillating  translating  :  -  v ,  and  the  general  case  of  the  radiation  (or  ship  motions) 
function  for  an  oscillating  translating  sour-  e.  cor  the  simple  ring  wave  case,  Liu  [1]  derived  near  and 
far  field  solutions  as  early  as  1952.  This  -  le  case  is  characterized  by  the  fact  that  there  is  exactly 
one  wave  at  each  field  point  and  in  the  far  field  these  waves  form  circular  rings  around  the  source. 
More  recently,  Grekas  and  Delhommeau  [2]  show  that  in  the  case  of  low  speed  the  far  field  ring 
waves  have  elliptic  shape.  The  far  field  wave  pattern  of  the  Kelvin  source  has  been  widely  studied 
since  this  gives  rise  to  wave  resistance,  the  drag  component  which  may  be  most  significantly  reduced 
by  suitable  ship  design.  It  is  well  known  that  there  are  two  waves  (transverse  and  divergent)  inside  a 
triangular  region  behind  the  ship,  with  a  half  angle  of  19.47°.  The  various  far  field  approaches  are 
summarized  in  the  classical  paper  by  Eggers,  Sharma,  and  Ward  [3].  It  is  only  in  recent  years  that 
Newman  [4],  and  Wang  and  Rogers  [5]  have  derived  efficient  computational  schemes  for  the  near 
field  flow. 

For  the  general  radiation  case,  where  the  source  is  both  oscillating  and  translating,  efficient 
robust  computational  schemes  for  the  near  field  flow  appear  to  be  lacking.  A  notable  exception  may 
be  the  work  by  Chang  [6].  In  the  far  field  case,  Newman  [7]  has  oudined  a  computational  approach 
for  obtaining  the  wave  pattern.  This  report  presents  a  numerical  implementation  of  this  formulation. 
First,  a  brief  description  of  the  derivation  of  the  far  field  expressions  is  given.  Next,  various  aspects 
of  carrying  out  their  numerical  evaluation  are  described  in  some  detaii.  These  include  an  efficient 
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iterative  scheme  for  finding  the  wave  propagation  directions  at  a  given  field  point,  the  treatment  of 
singular  points  and  lines  in  the  formulation,  a  grid  filtering  technique  to  eliminate  short  waves  which 
appear  as  noise,  and  an  efficient  method  of  extending  the  procedure  for  a  single  source  to  a  line  of 
sources  which  may  be  taken  to  represent  the  ship  hull  in  the  far  field.  Expressions  are  developed  for 
the  values  of  the  Green’s  function  and  its  slopes,  as  well  as  for  the  associated  wave  elevation  and  its 
slopes.  Numerical  results  are  presented  in  the  form  of  computer  generated  three-dimensional  and 
contour  plots  for  various  cases  of  a  single  source  and  a  line  of  sources.  These  figures  show  in 
particular,  the  differences  between  the  pattern  of  the  Green’s  function  and  the  wave  elevation, 
between  the  functions  (Green  and  elevation)  and  their  slopes,  and  between  the  wave  elevations  due  to 
a  single  source  and  a  line  of  sources  modeling  heave  (linear  vertical  oscillatory  motion)  and  pitch 
(rotational  oscillatory  motion  in  the  vertical  centerplane).  The  report  concludes  with  a  summary  of 
the  calculation  procedure  and  the  numerical  results. 

2.  DESCRIPTION  OF  FAR  FIELD  FORMULATION 

For  the  sake  of  uniformity,  the  present  report  will  consistently  follow  the  coordinate  system  and 
notation  used  in  [7].  This  coordinate  system  is  shown  in  Fig.  1,  and  corresponds  to  the  source  at  rest 
in  the  presence  of  a  current  c  in  the  —x  direction,  which  is  equivalent  to  the  source  traveling  at 
speed  c  in  the  +x  direction.  The  y  axis  corresponds  to  the  horizontal  direction  perpendicular  to  x , 
and  the  z  axis  is  directed  vertically  upwards. 

2.1  Derivation  of  Far  Field  Formulas  for  Values  of  Gu  for  Single  Souice 

Assume  that  the  Green’s  function  G  is  expressed  in  the  form 


G(x,y,z;£,77,f;r)  =  Re|cn(jr,y, fV"'  J 

where  x,y,z  and  are  respectively  the  coordinates  of  the  field  point  and  source  shown  in  Fig.  1, 
r  is  time,  and  w  is  the  circular  frequency  of  oscillation.  In  the  above  equation  and  in  the  remainder 
of  this  report,  only  the  real  part  of  the  complex  expression  is  of  interest  and,  for  the  sake  of  conveni¬ 
ence,  the  symbol  Re  will  be  omitted  hereafter.  Gu  must  satisfy  the  following  conditions 


V2Gn  =0;  z  <  0;  (x,y,z)  *  (£,*?, f) 

(2a) 

3Gu  ,  d2Gn  3Gu 

tu‘G  2iwc  +  c2  ,  +  g  .  ~  0 

dx  dx2  oz 

(2b) 

d  G„ 

Iim  =0 

z - oo  OZ 

(2c) 

lim  V  G  =  0 

(2d) 

R  —  06 


In  addition,  a  radiation  condition  of  no  incoming  waves  at  R  —  oo  is  needed  to  ensure  uniqueness  of 
the  solution.  Using  the  Fourier  transform  approach,  Gu  is  given  by 

Gn  =  7  ~  77  +  +  F2  +  F3  (3) 
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where  r  =  [(x  -  £)2  +  (v  “  0)2  +  (z  -  f)2]1'2,  r‘  =  [(x  -  $)2  +  (y  -  17)2  +  (z  +  D2]1/2. 

Ft,  Ft,  F3  are  double  integrals  of  the  complex  wavenumber  X  and  the  real  wave  propagation  direction 
u . 

Thf>  analysis  is  considerably  simplified  by  neglecting  terms  which  are  of  order  l/R  or  higher, 
where  R  is  related  to  x  and  y  by 


R  =  (x2  +  y2)1/2 

(4a) 

x  =  R  cos  0 

(4b) 

y  =  R  sin  0 

(4c) 

0  =  tan~l(y  /x) 

(4d) 

In  the  above,  0  is  the  spatial  direction.  By  using  the  calculus  of  residues  and  the  method  of  stationary 
phase  on  the  double  integrals  Fs,  F2,  and  F3,  the  following  algebraic  expression  for  the  far  field 
values  of  <7U  is  obtained 


Gu(R,d,z;£,ri ,f)  = 


8tt 

R 


1/2 


£  signfcos  u) 


X,  sin2  0 


sin2  u 


d0  /  m 
—  cos  ( u  -  6) 

du 


1/2 


m 


exp  "j  X, (u )[z  +  t  +  iR  cos  (u  —  0)  —  i(  cos  u  —  i rj  sin  u ]  +  —  +  sign 


^du2  j 


f}  (5> 


The  wavenumber  X,  (u )(/  =  1 ,2)  is  given  by 

X,  =  — ■  —  r —  [1  +  2r  cos  u  ±  (1  +  4r  cos  «)1/2]  (6) 

2  r  cos“  u 

where  the  upper  sign  is  used  for  i  =  1  and  the  lower  sign  for  /  =2 ,  t  ~  uc  /g ,  v  =  u>  2/g ,  g  is  the 
gravity  constant,  and  the  function  g,  (u )  is  given  by 


g,  =  X,  cos  (u  -  0)  (7) 

Equation  (5)  is  to  be  summed  over  the  values  of  u  which  give  dgj/du  =  0.  This  in  turn  gives  rise  to 
the  following  transcendental  u  -  0  equation 


ctn  0  =  - 


tan  u  ± 


(1  -f  4t  cos  u ) 
sin  u  cos  u 


(8) 
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where  u  is  in  the  interval  — x  <  u  <  0  -  x/2,  and  again  the  upper  and  lower  signs  are  to  be  used 
for  i  =  1  and  2,  respectively.  Fig.  2,  which  corresponds  to  Fig.  1  of  Ref.  [7],  shows  plots  of  the 
above  equation  for  various  values  of  the  parameter  r.  The  diagonal  line  connecting  the  points 
(0  =  0 ,  u  =  - x/2 )  and  (0  =  -*,u  =  x/2)  represents  the  upper  limit  of  permissible  values  of  u. 
This  figure  shows  that  there  may  be  anywhere  from  zero  to  5  roots,  depending  on  0  and  t. 


2.2  Limiting  Cases  of  r  —  0. 


Since  the  above  formulation  is  for  a  translating  oscillating  source,  it  should  contain  in  the  limit 
of  values  of  r(=ca>/g)  —  0  the  ring  wave  case  for  c  —  0  and  the  Kelvin  wave  case  for  «  —  0. 

The  ring  wave  case  is  recovered  by  taking  the  lower  sign  in  Eq.  (6),  setting  r  cos  u  =  e,  and 
expanding  the  square  root  term  to  second  power  in  e  resulting  in 


X2  - 


j"  1  +  2e  - 


(9) 


which  is  precisely  the  wavenumber  of  the  zero  speed  case.  The  Kelvin  waves  may  be  recovered  by 
taking  the  upper  sign  in  Eq.  (6),  noting  that  Wr2  =  u>2g2/c2u2  =  g2/c2  gives  the  fundamental 
wavenumber  \c  of  the  Kelvin  waves,  once  again  setting  r  cos  u  =  e,  and  expanding  the  square  root 
term  to  first  power  in  e,  resulting  in 


Xc 

2  cos2  u 


1  +  2e  + 

r 

1  +  T 

V  J 

. 

Xc  (2  +  2e)  Xc 
2  cos2  u  cos2  u  ’ 


0 


(10) 


The  term  Xc  /cos2  u  gives  the  precise  variation  of  the  wavelengths  of  the  Kelvin  wake  with  wave  pro¬ 
pagation  direction  u . 

2.3  Derivatives  of  the  Green’s  Function 

Often  it  is  of  interest  to  determine  not  only  the  values  of  the  Green’s  function  but  also  its  x  and 
y  derivatives.  In  the  case  of  damping  force  and  remote  detection  applications,  it  is  not  the  Green’s 
function  itself  but  the  associated  wave  elevation  and  slopes  which  are  of  principal  interest.  All  of 
these  derivative  quantities  may  be  conveniently  obtained  from  the  formula  for  Cu  given  in  Eq.  (5). 

Recalling  that  this  equation  is  a  far  field  expression  for  large  values  of  R,  the  operations 
d/dx  =  d/d(R  cos  0)  and  d/dy  =  d/d(R  sin  0)  would  add  a  power  of  l/R  to  most  terms  in  Eq.  (5). 
The  only  term  for  which  this  is  not  the  case  involves  the  exponential  function 


r 

d  d 

d  d 

to'  dy  J 

exp  [X,  i7?  cos  (u  -  0)]  = 

dx'  dy  ^ 

[X,  i  ( x  cos  u  +  y  sin  u )] 


(U) 


for  which  the  derivatives  d/dx  and  d/dy  are  equivalent  to  multiplying  the  function  by  iX,  cos  u  and 
iXj  sin  u ,  respectively.  Thus,  the  x  and  y  derivatives  of  Gu  are  simply  given  by 

1/2 

£  iXj  cos  u  G||,  (12a) 


dGn  =  8x 
dx  R 


4 


3£ n 
dy 


8tt 

/? 


1/2 

£  lX< 


sin  «  G11( 


(12b) 


where  GUj  refers  to  the  entire  expression  after  the  summation  sign  in  Eq.  (5).  The  elevation  zs  of 
the  wave  pattern  for  the  present  case  of  an  oscillating  source  in  the  presence  of  a  current  c  in  the  —x 
direction,  shown  in  Fig.  1,  is  given  by 


-1 

^Gu 

<?GU] 

f  8tt 

*  *  ! 

dt 

dx  J 

L*. 

1/2 


£  -  -  c\  cos  u)G\\i 

g 


(13) 


The  x  and  y  slopes  are  then  obtained  similar  to  Eqs.  (12) 


1 

r?  | 
<X5 

dx  1 

L*  J 

1 

1 

j 

r 

87 r 

dy 

R 

J 

1/2 


V,  — -  i X.  cos  u(u  -  c\  cos  u)  GUl 
g 


1/2 


-l  X 


V  — -  iX.  sin  u(oi  —  cX.  cos  u )  GUl 
8 


(14a) 


(14b) 


2.4  Extension  to  Multiple  Sources 

The  straightforward  way  of  evaluating  the  Green’s  function  for  a  series  of  N  sources  with 
strengths  Q;  and  locations  (£;  ,i?7 ,  )  is  to  repeat  the  calculation  of  Eq.  (5)  N  times  and  sum  over  the 

results.  Actually,  as  a  result  of  the  far  field  assumptions  inherent  in  Eq.  (5),  the  computational  effort 
for  N  sources  differs  very  little  from  that  for  one  source.  Most  of  the  computational  effort  in 
evaluating  Gu  occurs  in  the  iterative  search  for  the  various  u -roots  of  Eq.  (8)  for  a  given  value  of  8. 
Eq.  (5)  shows,  however,  that  8  and  R  are  always  measured  from  the  origin  even  though  the  source 
may  be  located  at  ($,17, f).  This  assumption  is  possible  since  R  »  (£,17, f).  Thus,  only  the  exponen¬ 
tial  term  exp(X,(u)[f  -  i  £  cos  u  -  it;  sin  «]j  varies  for  the  different  sources.  The  term  X,  f  gives 
the  decrease  in  the  strength  of  the  source  due  to  submergence  f  below  the  free  surface  while  the  term 
-i£  cos  u  -  irf  sin  u  gives  the  shift  in  the  phase  of  the  wave  originating  from  (£,77)  relative  to  the 
wave  due  to  a  source  located  at  (0,0). 


A  particularly  useful  configuration  for  modeling  ship  hulls  in  the  far  field  is  a  line  of  N  sources 
with  the  same  values  of  17  and  f,  but  different  values  of  Qj  and  £;  .  In  this  case,  the  value  of  Gu  at  a 
given  field  point  ( x  ,y  ,z )  is  given  by  the  following  double  sum 


G„  = 


87T 

R 


1/2 


Z  sign  (cos  u)  J 


\  sin2  8 


•  2  1  1 

sin  u  |  ——  cos  (u  -  8)  \ 

du 


1/2 


X 


exp 


+  f  +  *7?  cos  (u  -  6)  -  i  tj  sin  u  ]  +  —  +  sign 


x 


5 


(15) 


v  a 

Y.  —  exp  (-i£.  cos  u) 

,  =  ,  4x; 


Similar  expressions  hold  for  extending  the  derivatives  of  Cu,  the  elevation  zs  and  its  slopes  for  a  sin¬ 
gle  source,  given  in  Eqs.  (12)  to  (14),  to  a  line  of  sources. 

3.  COMPUTER  IMPLEMENTATION  OF  THEORETICAL  FORMULATION 

The  preceding  formulation  has  been  numerically  implemented  in  the  form  of  the  FORTRAN  77 
computer  program  RADGREEN.  The  present  chapter  describes  the  more  complex  or  novel  aspects 
of  the  implementation.  Input,  output,  and  computer  time  requirements  of  the  program  are  also 
described  to  give  an  indication  of  actual  usage. 


3.1  Iterative  Solution  for  Wave  Propagation  Direction 

The  most  time  consuming  part  of  the  procedure  is  the  finding  of  the  zero  to  five  values  of  u 
which  correspond  to  a  given  value  of  6  in  Eq.  (8).  The  solution  scheme  essentially  consists  of  a  glo¬ 
bal  search  to  bracket  the  root  values  of  u  within  certain  bounds,  and  then  a  localized  search  to  pin¬ 
point  the  precise  values. 

The  global  scheme  is  contained  in  subroutine  FINROOT,  which  evaluates  the  error  functions  /  j 
and  f2  (for  given  values  of  6  and  r)  for  721  equally  spaced  values  of  u  between  - 180  and  0  degrees. 
The  error  functions  / ,  and  f2  are  defined  as  the  difference  between  the  left  and  right  sides  of  Eq.  (8) 
using  the  upper  (X,  case)  and  lower  (X2  case)  signs,  respectively.  Where  the  function  changes  sign 
between  two  adjacent  values  of  u ,  indicating  a  zero  crossing  and  hence  the  location  of  a  root,  these 
values  of  u  are  identified  and  saved  as  initial  bounds  for  a  more  precise  determination  by  subroutine 
FR1.  This  subroutine  closes  in  on  the  value  of  the  root  by  an  interval  halving  technique  where  the 
error  function  is  evaluated  at  the  midpoint  between  the  bounds,  and  then  moving  the  bound  with  the 
same  sign  of  the  error  function  to  the  midpoint.  This  iterative  procedure  is  continued  until  the  error 
function  has  an  absolute  value  less  than  10-6.  The  resulting  value  of  u  and  the  associated 
wavenumber  case  (X)  or  X2)  are  then  returned  to  the  main  program  for  use  in  Eq.  (5). 

3.2  Treatment  of  Singular  Points  and  Lines 

Perhaps  the  most  obvious  indication  that  Eq.  (5)  is  for  the  far  field  is  the  factor  1  r^R  ,  which 
—  oo  as  R  —  0.  Where  calculations  are  desired  for  field  points  relatively  close  to  the  origin,  the  pro¬ 
gram  provides  the  option  for  the  user  to  input  a  minimum  bound  on  values  of  R  to  be  used  in  Eq. 
(5).  For  the  zero  speed  case,  Liu  [1]  shows  that  the  far  field  formulation  becomes  valid  for 
R  >  0.3  /2,  where  l2  is  the  wavelength  and  is  related  to  the  wavenumber  X2  given  in  Eq.  (9)  by 


2x  _  2 icg 

X2  op- 


(16) 


Eq.  (5)  shows  that  the  factor  contained  within  brackets  may  be  singular  at  values  of  8  where 
sin  u  =  0,  cos  («  —  6)  =  0,  dd/du  =  0,  and  (due  to  Eq.  (6)  for  X()  cos  u  =  0.  Actually,  the 
singularities  due  to  sin  u  and  cos  u  are  completely  removable  while  the  behavior  at  cos  (u  -  8)  —  0 
is  also  regular.  Thus,  only  the  case  dd/du  =  0  is  an  actual  singularity  in  the  calculation. 
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The  removable  nature  of  the  singularities  for  sin  u  and  cos  u  may  best  be  illustrated  by  explicity 
writing  out  dd/du  obtained  by  differentiating  Eq.  (8)  with  respect  to  u ,  resulting  in 


d6_ 

du 


sin:  6 


sin'  u  cos  u 


- — [sin'  u  ±  (  — 2r  sin  u )  sin  u  cos  u  1  +  4 r  cos  u 


±(-l)Vl  +  47  cos  u  cos  2m] 


(17) 


It  is  seen  that  the  factors  sin2  6/ sin2  u  and  1/cos2  n  are  cancelled  by  corresponding  factors  in  Eqs. 
(5)  and  (6).  Thus,  Eq.  (5)  is  entirely  free  of  singularities  in  terms  of  possible  zeroes  of  sin  u  and 
cos  u . 

The  regular  behavior  at  cos  (u  -  6)  —  0  may  be  seen  by  noting  that  this  implies  (since 
|  u  —  6 1  =  t/2  or  3ir/2  in  Eq.  (8))  that 


1  +  4t  cos  u,  =  0 


(18) 


and,  from  Eq.  (6),  that  \!  =  X2.  From  Fig.  2,  and  also  by  considering  Eq.  (17),  it  is  seen  that  at 
these  values  of  u  the  slope  du/dd  =  0,  so  that  the  product  cos (m  -  6)d9/du  has  the  indeterminate 
behavior  0/0.  Numerical  calculations  show  thai  the  product  has  finite  values,  with  resultant  values  of 
the  magnification  factor  n  decreasing  from  7.0  for  r  =  0.26  to  0.5  for  r  =  1.0,  where  n  is  given  by 


H  - 


sin2  6 


sin2  m 


d6_ 

du 


cos  (m  -  6) 


(19) 


The  calculations  show  that  the  values  of  n  vary  relatively  slowly  with  small  changes  in  u  away  from 
u, .  Since  the  indeterminate  behavior  does  not  involve  two  factors  which  exactly  cancel  as  in  the  pre¬ 
vious  cases  of  sin  u  and  cos  u ,  it  is  numerically  necessary  to  give  special  treatment  to  this  factor  for 
u  near  u, .  This  is  implemented  in  the  computer  program  by  restricting  u ,  for  cases  of  t  >  0.25  and 
m  <  0,  to  be  at  least  0.1  deg  greater  than  u, ,  i.e.  u  >  u,  4-0.1  deg.  To  obtain  the  correct  limiting 
behavior,  the  values  of  6  corresponding  to  this  revised  value  of  u  must  be  obtained  from  Eq.  (8)  and 
used  in  Eq.  (5). 

Figure  2  shows  that  for  6  near  180  deg,  there  may  be  two  values  of  6,  8C,  at  which  there  is  a 
vertical  tangent  dd/du  =  0.  These  are  the  so-called  cusp  lines  which  give  rise  to  unbounded  values 
of  G,,  and  also  serve  as  bounding  lines  such  that  the  waves  are  contained  within  the  angular  region 
6C  <  6  <  180  deg.  This  singular  behavior  illustrates  the  limitations  of  the  present  method  since 
actual  values  of  Gu  are  finite  at  the  cusp  lines  and  also  extend  outside  these  lines.  Higher  order 
asymptotic  methods,  such  as  that  of  steepest  descent,  are  currently  under  investigation  to  obtain  more 
accurate  calculations  of  Gn  on  and  near  the  cusp  lines. 

A  particularly  simple  and  well  known  case  is  obtained  by  setting  Eq.  (17)  equal  to  0,  with  r  = 
0,  and  taking  the  upper  sign  (Xt  case)  resulting  in  the  following  equation  for  uc 

tan2  uc  =  1/2  o  uc  =  35.16  deg  (20) 
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By  substitiKij>r!  mis  equation  in  Eq.  (8)  the  corresponding  value  of  0C  is  found  to  be 

6C  =  180.  -  19.47  =  160.53  deg  (21) 

This,  of  course,  refers  to  the  well  known  Kelvin  waves  for  a  nonoscillating  translating  ship  for  which 
the  transverse  and  divergent  waves  meet  at  the  cusp  line  with  propagation  angle  35.16  deg  and  the 
waves  are  confined  within  an  angle  at  19.47°  to  the  —x  axis.  Since  u  >  0,  these  waves  (located 
above  the  u  =  0  axis  in  Fig.  2)  are  propagating  against  the  current  (i.e.,  with  the  ship)  while  the 
remaining  waves  (located  below  the  u  =  0  axis)  are  propagating  with  the  current. 

It  is  of  interest  to  ascertain  the  manner  of  this  singular  behavior  as  u  —  uc  for  /  =  0.  The 
close  bunching  of  the  different  curves  for  the  Kelvin  like  waves  above  the  u  =  0  axis  suggests  that 
this  behavior  may  also  be  applicable  for  other  values  of  r.  The  analysis  consists  of  evaluating  the 
magnification  factor  /t  for  small  increments  A u  away  from  uc ,  i.e.,  evaluating  Eq  (19)  with 
u  =  uc  +  A u  .  Numerical  evaluation  of  the  subsequent  expression  shows  that  n  varies  approximately 
as 


H  =  13.5/Aw  (22) 

where  A u  is  expressed  in  degrees.  For  a  number  of  reasons  it  is  necessary  to  set  an  upper  limit  on  /x. 
Besides  possible  numerical  overflow  problems,  excessively  large  values  near  the  cusp  line  would  tend 
to  overwhelm  the  values  elsewhere  so  that  they  would  appear  to  be  nearly  zero  in  plots  of  the  func¬ 
tion.  In  the  computer  program  the  default  value  of  n  is  set  equal  to  100.  This  corresponds  to  taking 
A u  =  0.135  deg  in  Eq.  (22).  Taking  a  lower  value  of  /i  =  10,  to  give  more  prominence  to  off-cusp 
values  of  Gn,  corresponds  to  a  value  A u  -  1.35  deg. 

3.3  Grid-Based  Filtering  of  Short-Wavelength  Noise 

Figure  2  shows  that  the  wavelengths  change  with  0  and,  in  general,  a  wide  range  of  wavelengths 
may  be  expected  for  an  arbitrary  x  —  y  grid  of  field  points.  Those  waves  where  lengths  are  shorter 
than  or  comparable  to  the  lengths  Ax  or  Ay  of  the  individual  grid  cell  will  show  up  as  high  frequency 
noise.  In  order  to  remove  this  noisy  behavior  in  plots  of  the  wave  field,  the  program  omits  calcula¬ 
tions  for  wavenumbers  X,(/  =  1,2)  which  exceed  the  maximum  wavenumber  \M  (minimum 
wavelength  lm )  given  by 


n  [min(Ax  ,Av)J 

where  n  is  the  number  of  grid  points  used  to  model  the  minimum  wavelength.  The  default  value  is 
4.  A  higher  value  of  n  may  be  used  to  focus  on  only  the  longer  waves,  while  a  lower  value  of  n 
may  be  used  to  account  for  the  shorter  waves  in  spite  of  their  possible  noise-like  appearance. 

3.4  Description  of  Computer  Program 

Program  RADGREEN  consists  of  a  main  program  and  the  subroutines  CPTIME,  FINROOT, 
FRI,  and  UNFORM.  The  main  program  accepts  input  data,  performs  various  initializations,  and  car¬ 
ries  out  the  calculations  given  in  Eqs.  (4)-(15)  for  a  rectangular  x  -  v  grid  of  field  points.  The 
CPTIME  subroutine  keeps  track  of  the  actual  CPU  time  used  for  the  calculations.  As  mentioned  pre¬ 
viously,  the  subroutines  FINROOT  and  FRI  respectively  perform  global  and  local  searches  for  the 
roots  u  of  Eq.  (8).  The  subroutine  UNFORM  writes  onto  an  unformatted  output  file  the  calculated 
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aiue  of  the  Green’s  function  Gu  (or  the  wave  elevation  and  its  slopes  in  a  form  suitable  for 
three-dimensional  and  contour  plots. 


Input  data  are  entered  by  means  of  an  input  file  which  conforms  to  the  following  READ  state¬ 
ments. 

READ  (5,521)  filename 
521  FORMAT  (A30) 

READ  (5,*)  '•’I,  gcal,  nom,  nsr,  ip:,  wnp,  rmf,  fmumx 
READ  (5,*)  (omga  (m),  m  =  1,  nom) 

READ  (5,*)  lcol,  krow,  xO,  xdel,  yO,  ydel,  eta,  zeta 

DO  999  iom  =  1 ,  nom 

READ  (5,*)  (qamp  (m),  m  =  1,  nsr) 

READ  (5,*)  (qphd  (m),  m  =  1,  nsr) 

READ  (5,*)  (qxi(m),  m  =  1,  nsr) 

The  input  variables  are  defined  as  follows: 

1.  filename  is  the  output  file  onto  which  are  written  the  calculated  values. 

2.  vel  is  the  current  or  ship  speed  c  in  ft/s. 

3.  gcal  is  a  calculation  indicator  such  that  for  gcal  >0.1,  Gn  is  calculated;  otherwise,  zs  is 
calculated. 

4.  nom  is  the  number  of  values  of  u  to  be  considered,  1  <  nom  <11. 


5.  nsr  is  the  number  of  sources  on  the  line  1  S  nsr  <  21. 

6.  ipr  is  an  output  indicator  such  that  for  ipr  >  1  the  calculated  function  (Gu)  or  zs  and  its 
slopes  are  written  onto  the  output  file;  otherwise,  only  the  function  itself  is  written. 

7.  wnp  is  n ,  defined  in  Eq.  (23),  the  number  of  grid  points  used  to  model  the  minimum 
wavelength  which  will  be  considered  in  the  calculations;  entering  wnp  <  0.  triggers  the 
default  value  wnp  =  4. 

8.  rmf  is  the  multiplier  to  establish  Rmin  (the  minimum  value  of  R)  to  be  used  in  Eq.  (5) 
where  R ^  =  rmf  *  min  (Ax, Ay);  if  rmf  is  input  <  0.,  the  default  value  rmf  =  10.  is 
used. 

9.  fmumx  is  the  maximum  value  of  the  magnification  factor  fi,  defined  in  Eq.  (19);  entering 
fmumx  <  0.  triggers  the  default  value  fmumx  =  100. 

10.  omga  (m)  is  the  mth  value  of  u>  in  rad/s,  1  <  m  <  nom. 

11.  lcol  is  the  number  of  grid  points  in  the  x  direction,  1  <  lcol  <  151. 

12.  krow  is  the  number  of  grid  points  in  the  y  direction,  1  <  krow  <  151. 

13.  xO  is  the  value  of  x  in  ft  for  the  initial  grid  point. 
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14.  xdel  is  the  value  of  Ax  in  ft,  the  size  of  each  grid  cell  in  the  x  direction. 

15.  vO  is  the  value  of  y  in  ft  for  the  initial  grid  point. 

16.  ydel  is  the  value  of  Ay  in  ft,  the  size  of  each  grid  cell  in  the  y  direction. 

17.  eta  is  j?  in  ft,  the  y  value  of  the  location  of  the  line  of  sources. 

18.  zeia  is  f  in  ft,  the  z  value  of  the  location  of  the  line  of  sources. 

19.,  20.,  21.  qamp  (m).  qphd  (m),  and  qxi  (m)  are  respectively  the  amplitude  in  ft3/s,  the  phase 
in  deg,  and  the  x  -location  in  ft  of  the  mth  source,  1  <  m  <  nsr;  repeat  for  each 
new  value  of  omga. 

It  should  be  emphasized  that  the  singularity  moderating  factors  rmf  and  fmumx,  and  the  noise 
suppression  factor  wnp,  have  been  implemented  mainly  to  enhance  the  smoothness  and  visibility  of 
overall  plots  of  the  calculated  functional  values.  If  such  considerations  are  not  of  concern  to  the  user, 
and  the  actual  values  of  Eq.  (5)  are  desired,  then  simply  enter  arbitrarily  small  nonzero  values  for  rmf 
and  wnp  (such  as  0.001)  and  a  large  value  for  fmumx  (such  as  9999). 

Computer  time  at  one  field  point  for  a  single  value  of  w  depends  largely  on  the  given  value  of  r. 
Using  the  Hewlett  Packard  9000,  Model  550  minicomputer,  CPU  time  typically  decreases  from  0.3 
seconds  for  r  ~  0  to  0.2  seconds  for  r  =  1 .  Corresponding  total  CPU  time  for  the  maximum  grid 
of  151  x  151  field  points  would  vary  between  115  and  75  minutes,  approximately.  Due  to  the 
approach  used  in  Eq.  (15),  the  calculation  time  for  a  line  of  sources  differs  little  from  that  for  a  single 
source. 

4.  NUMERICAL  RESULTS 

A  number  of  calculated  results  are  presented  to  illustrate  the  capability  of  the  program  and  the 
types  of  wave  patterns  to  be  expected  from  a  single  source  and  a  line  of  sources  for  different  values 
of  the  fundamental  parameter  t  =  cui/g .  These  results  are  presented  in  the  form  of  three-dimensional 
and  contour  plots.  Most  of  the  results  are  presented  in  the  form  of  the  wave  elevation  z,  given  by 

Eq.  (13).  Some  results  are  also  presented  for  the  Green’s  function  itself,  Eq.  (5),  as  well  as  the 

slopes  defined  by  Eqs.  (12)  and  (14). 

By  way  of  illustrating  the  symmetry  properties  about  the  x-axis,  Figs.  3a,  b,  and  c  respectively 
sh  sample  plots  of  the  wave  elevation  zs ,  the  x  slope,  and  y  slope  for  t  =  0.125  for  the  square 

-500  <  x  <  500  and  -500  <  y  <  500.  These  figures  show  that  the  calculated  results  con- 

the  expected  symmetry  of  zs  and  its  x  slope,  and  the  antisymmetry  of  the  y  slope  about  y  =  0. 

which-  ,“>ay  be  deduced  from  Eqs.  (5),  (13),  and  (14).  Similar  symmetry  properties  hold  for  the 
Gr te..  -  function  and  its  slopes. 

In  view  of  the  above  symmetry  properties  about  y  =  0,  the  remaining  figures  are  shown  only 
for  one  side  of  y,  0  <  y  s  450  ft,  and  -450  s  x  <  450  ft,  with  Ax  =  Ay  =  6  ft.  The  line  of 
sources  is  taken  to  be  submerged  6  ft  (f  =  -  6)  and  with  y  coordinate  =  0  (v  =  0).  Calculations 
are  performed  for  a  single  source  with  x  coordinate  =  0  ({  =  0)  and  for  10  equally  spaced  sources 
with  -135  <  £  S  135  ft.  This  line  may  be  taken  to  be  the  far  field  representation  of  a  ship  with 
draft  12  ft  and  length  300  ft.  The  strengths  of  the  10  sources  are  assumed  to  be  either  (a)  equal  or 
(b)  linearly  proportional  to  the  f  location  of  the  source.  Case  (a)  is  an  approximation  to  heave 
motions  where  the  ship  oscillates  in  the  z  direction,  while  case  (b)  is  an  approximation  to  pitch 
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motions  where  the  ship  undergoes  rotational  oscillations  about  the  y  axis.  The  angular  frequency  a>  is 
fixed  at  0.805  rad/s,  which  is  representative  of  Sea  State  5  conditions. 

To  show  the  difference  between  the  wave  elevation  zs  and  Green’s  function  Gn  cases.  Figs.  4a, 
b,  and  c  respectively  show  contour  plots  of  zs,  dzs/dx,  and  dz5/dy  for  r  =  0.125,  while  Figs.  5a  to  c 
show  corresponding  plots  for  Gn.  These  figures  show  that  the  differences  are  relatively  small. 

Figures  6a  to  f  show  three-dimensional  plots  of  zs  for  the  heave  case  for  ship  velocity  c  =  0.2, 
5.0,  9.8,  10.2,  20.0,  and  40.0  ft/s  which  correspond  to  r  =  0.005,  0.125,  0.245,  0.255,  0.5,  and 
1.0.  The  highest  speed  may  be  considered  representative  of  the  design  speed.  These  figures  show 
how  the  wave  pattern  changes  from  ring  waves  to  Kelvin  type  waves  behind  the  ship  as  the  velocity 
increases.  Also,  the  previously  mentioned  cusp  lines,  which  correspond  to  values  of  6  giving  the 
vertical  tangents  in  Fig.  2,  appear  for  values  of  t  >  0.245. 

To  illustrate  the  effect  of  source  distribution  on  the  wave  pattern.  Figs.  7a,  b,  c  respectively 
show  contour  plots  of  zs  for  c  =  0.2  (r  =  0.005)  for  the  single  source,  heave,  and  pitch  motions, 
while  Figs.  8a  to  c  show  corresponding  plots  for  c  =  9.8  (t  =  0.245).  Figure  7a  shows  the 
expected  symmetry  about  x  =  0  of  the  ring  waves  from  a  single  source  at  this  near  zero  value  of  c 
Relative  to  the  nearly  circular  waves  of  the  single  source,  Fig.  7b  shows  that  the  heave  motion  tends 
to  concentrate  the  waves  to  the  side  of  the  ship  (along  the  y  axis)  while  the  pitch  motion  tends  to 
gives  the  largest  waves  fore  and  aft  of  the  ship  (along  the  x  axis).  In  the  case  of  the  the  more  com¬ 
plex  wave  pattern  for  c  =  9.8  (r  =  0.245),  the  single  source  and  pitch  cases  give  similar  wave  pat¬ 
terns  which  are,  in  turn,  more  varied  than  those  for  the  heave  case. 

5.  SUMMARY 

A  previously  formulated  approach  for  calculating  the  far  field  behavior  of  the  radiation  Green’s 
function  Gu,  corresponding  to  a  source  which  is  translating  and  oscillating  beneath  the  free  surface, 
has  been  numerically  implemented  in  the  form  of  computer  program  RADGREEN.  Supplementary 
expressions  have  been  derived  for  the  far  field  derivatives  of  Gn  as  well  as  for  the  associated  wave 
elevation  and  is  derivatives.  The  program  first  uses  a  global  scheme  and  then  a  localized  scheme  to 
determine  the  wave  propagation  directions  u  for  a  given  spatial  angle  0  of  the  field  point.  It  is  shown 
that  the  behavior  is  regular  at  several  apparent  singular  values  of  0  and  u.  However,  at  cusp  lines 
which  correspond  to  locations  of  vertical  tangents  in  Fig.  2,  the  function  does  become  unbounded  and 
it  is  necessary  to  impose  limits  on  the  magnification  of  the  waves  at  these  points.  It  is  also  shown 
that  the  grid  filtering  technique  removes  short  wavelength  noise  by  retaining  only  those  waves  which 
have  lengths  equal  to  several  cells  of  the  computation  grid.  Calculation  time  on  the  HP  9000,  Model 
550  minicomputer  for  one  field  point  varies  from  0.3  seconds  to  0.2  seconds,  depending  on  the  value 
of  the  dimensionless  parameter  t  =  cu/g . 

A  series  of  three-dimensional  and  contour  plots  are  given  to  illustrate  the  capability  of  the  pro¬ 
gram  and  to  give  representative  wave  patterns  for  a  single  source  and  for  a  ship  undergoing  heave  and 
pitch  motions.  The  heave  and  pitch  motions  are  obtained  by  approximating  the  ship  hull  by  a  series 
of  10  sources  on  a  line  with  source  strengths  which  are  respectively  equal  to  each  other  and  propor¬ 
tional  to  the  longitudinal  location  of  the  source.  These  plots  show  that  the  program  reproduces 
expected  symmetry  properties  about  the  x  and  y  axes.  There  is  relatively  little  difference  between 
contour  plots  of  Gn  and  its  derivatives  and  corresponding  plots  of  zs  and  its  derivatives.  At  low 
speeds,  relative  to  a  single  source,  the  wave  pattern  for  heave  motion  tends  to  concentrate  the  waves 
to  the  side  of  the  ship  while  the  pitch  motion  produces  the  largest  waves  fore  and  aft  of  the  ship. 
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FREE  SURFACE 


•  (x,  y,  z)  Field  Point 


CURRENT  c 


•  i£fV> ft  Source  Point 


Fig.  1  —  Definition  of  coordinate  system 
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Curves  representing  the  u  -  6  Equation  (8)  for  various  values  of  r 


Fig.  3a  —  Elevation  zs 


Fig.  3b  —  ;c -slope  dzs/dx 

Fig.  3  —  Sample  three-dimensional  plots  of  wave  patterns,  r  =  0. 125 
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y 


Fig.  3c  —  y -slope  bzs  /by 

Fig.  3  —  (Continued)  Sample  three-dimensional  plots  of  wave  patterns,  r  =  0.125 


Fig.  4a  —  Elevation  zs 


"ig.  4  —  Contour  plots  of  wave  elevation  zs  for  single  souce,  t  =  0.125 
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Fig.  5a  —  Green’s  function  Gn 


Fig.  5b  —  * -derivative  dGu/dx 

Fig.  5  —  Contour  plots  of  Green’s  function  Gn  for  single  source,  r  =  0.125 
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Fig.  5c  —  ^-derivative  dGu/dy 

Fig.  5  —  (Continued)  Contour  plots  of  Green’s  function  Gu  for  single  source,  r  =  0.125 


x 


Fig.  6a  —  c  =  0.2  ft/s,  r  =  0.005 
Fig.  6  —  Three-dimensional  plots  of  heave  wave  elevation  zs 
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Fig.  6b  —  c  =  5.0  ft/s,  r  —  0.125 


Fig.  6c  —  c  =  9.8  ft/s,  r  =  0.245 


Fig.  6  —  (Continued)  Three-dimensional  plots  of  heave  wave  elevation  zs 
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Fig.  6f  —  c  =  40.0  ft/s,  r  =  1.0 

Fig.  6  —  (Continued)  Three-dimensional  plots  of  heave  wave  elevation  zs 
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Fig.  7a  —  Single  source 


Fig.  7b  —  Heave  for  line  of  sources 
Fig.  7  —  Contour  plots  of  wave  elevation  zs,  c  =  0.2  ft/s,  r  =  0.005 
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Fig.  8a  —  Single  source 

Fig.  8  —  Contour  plots  of  wave  elevation  zs,  c  =  9.8  ft/s,  r  =  0.245 
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Fig.  8b  —  Heave  for  line  of  sources 


Fig.  8c  —  Pitch  for  line  of  sources 

Fig.  8  —  (Continued)  Contour  plots  of  wave  elevation  zs ,  c  =9.8  ft/s,  r  =  0.245 


